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Abstract. Pleistocene climate oscillations had a dramatic impact on the distribution of temperate ectothermic organisms. 
These range fluctuations plausibly have left a footprint on the species’ genetic structure which is linked with modelled 
paleoclimatic conditions, allowing us to infer which environmental factors shaped the evolutionary history of Sand liz- 
ards. In this study, we evaluated the patterns of niche diversification of Lacerta agilis in Eurasia, based on mitochondrial 
DNA analyses and ecological niche models. The lineage of Sand lizards evolved 3.8-1.8 Mya, being the most basal in the 
east-Caucasian subclade, with eight mitochondrial subclades separated into two groups. These groups represent two in- 
dependent waves of expansion from the ancestral region constituted by the Caucasus and the adjacent northern plains to 
the Paleartic. Paleoclimatic models suggested a high instability of the range of this lineage in the last 3 Mya, with niche 
contractions during the colder glacial periods and expansions following the glacier retreat. This suggested an allopatric di- 
versification process, with subspecies boundaries upon secondary contact maintained by competitive interactions, at least 
between closely related pairs. The great mountain systems of the Mediterranean Peninsulas constituted stable refugia dur- 
ing the Pliocene-Pleistocene cycles, favouring the evolution of endemic subclades. These montane subclades have higher 
mitochondrial diversity than those that occur in the plains. However, the Pyrenean endemic L. agilis garzoni is exception 


possibly due to the occupation of a very small refugium during the recent glacial phases. 
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Introduction 


Recent glacial and interglacial cycles have played a crucial 
role in biotic diversification in the temperate regions of the 
northern hemisphere (HEWITT 2000). Ectothermic organ- 
isms such as reptiles are suitable models to assess the ef- 
fects of Pleistocene climatic oscillations because of their 
entire dependence on external sources to thermoregulate 
and lower ability to disperse compared to mammals and 
birds (PAULO et al. 2001, 2008, URSENBACHER et al. 2008). 
Recent studies locate centres of higher lineage diversifica- 
tion of Palearctic reptiles in the great mountain ranges of 
southern Europe, the Caucasus, Carpathians, Dinaric Alps, 
Pindus, Apennines, and Pyrenees (JOGER et al. 2007). This 
suggests that these mountain ranges were stable climate 
refugia through the Pleistocene. 

The Sand lizard (Lacerta agilis) is one of the more wide- 
ly distributed reptiles of Eurasia, occurring from Mongo- 
lia through Central Asia to the western European Atlan- 
tic plains, including the British Isles (BIscHOFE 1988). This 


species could have an intricate biogeographic history, as in- 
dicated by previous phylogenetic studies (KALYABINA et al. 
2001, ANDRES 2014). This genetic structure is also associ- 
ated with some differences in the external morphology that 
has led to the description of several subspecies based on 
the geographical variation of these traits (BISCHOFF 1988). 

One of these subspecies is L. agilis garzoni, a micro en- 
demism confined to the eastern Pyrenees (AMAT et al. 
2003). Lacerta a. garzoni is distributed allopatrically with 
L. agilis agilis, which is widespread across the western Eu- 
ropean Atlantic Plain, reaching north to south of the Scan- 
dinavian Peninsula (BISCHOFF 1988). This latter subspecies 
intergrades across a narrow zone with the L. a. argus along 
central Europe, which is replaced in the eastern margins 
of the European Atlantic Plain by L. a. chersonensis. The 
mountains of the Balkans are inhabited by another en- 
demic subspecies, L. a. bosnica. Lacerta a. exigua contacts 
L. a. chersonensis along the plains of Ukraine, Bielorussia, 
and Russia, extending further to the east, central Asia, and 
reaching Mongolia. The mountains around the Black Sea 
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are occupied by other endemic subspecies: L. agilis tauridi- 
ca (Crimea, Ukraine) and L. a. boemica (eastern Caucasus, 
KUKUSHKIN et al. 2020). 

Based on morphological characters and molecular data, 
previous studies indicated that this species originated on the 
Caucasus-Black Sea coast during the late Miocene to early 
Pliocene (YABLOKOV 1981, KALYABINA et al. 2001, JOGER et 
al. 2007). The current species range and the complex lineage 
genetic structure and morphological variability could be in- 
fluenced by tectonics and paleoclimate conditions (KALYA- 
BINA et al. 2001, JOGER et al. 2007). These studies suggested 
an early split between western and eastern basal lineages 
during the late Pliocene, followed by several waves of ex- 
pansion from the Caucasus region to western Europe. The 
Balkans and Carpathian mountains were the glacial refugia 
from which the species reached central Europe not earlier 
than 100,000 years with a more recent (Holocene) coloni- 
zation of the Pyrenees (KALYABINA et al. 2001). 

This complex diversification process could be governed 
by niche conservatism, as indicated by the relict presence of 
several montane subclades in the Mediterranean peninsulas. 
Niche conservatism has been invoked to explain patterns of 
speciation and niche diversification in several lineages of 
temperate reptiles, possibly driven by the Pliocene-Pleis- 
tocene climatic instability (SAGoNAs et al. 2014, JABLONSKI 
et al. 2019). According to this hypothesis, the species would 
have expanded geographically along favourable environ- 
mental corridors, whose extension fluctuated following the 
alternation of warm-cold periods. In this scenario, sepa- 
rated subclades only occupy climatically analogous regions 
throughout the species’ vast area of distribution. The range 
boundaries between parapatric populations do not follow 
any environmental discontinuity and are better explained 
by competitive interactions (FREEMAN 2015). A second test- 
able hypothesis is that this vast geographic range could be 
explained by niche shifts. In this scenario, some subclades, 
perhaps those geographically distant or those with more 
deep genetic divergence, have adapted to use new environ- 
ments (PYRON & BURBRINK 2009). In this case, the separa- 
tion between parapatric taxa or lineages clearly follows an 
environmental discontinuity. In addition, this process of 
niche diversification would generate a weaker phylogenetic 
signal in the environmental niche than the previous one. 

The goal of our research is to evaluate the following ques- 
tions: a) determine whether the L. agilis mitochondrial ge- 
netic structure is associated with past environmental cycles; 
b) assess the possible diversification mechanisms of this lin- 
eage: conservatism versus niche shifts and which variables 
are shaping niche diversification; c) describe the evolution of 
the niche of L. a. garzoni, including the environmental fac- 
tors which may have determined its micro-endemic status. 


Material and methods 
Phylogenetic inference and biogeography 


We used DNA samples of 35 individuals from nine locali- 
ties along the Pyrenean range of L. agilis and three samples 


from Germany, Hungary, and the United Kingdom stored 
in the Institut de Biologia Evolutiva (CSIC-UPF) (Sup- 
plementary document S1). Molecular sampling was per- 
formed using the mitochondrial gene cytochrome b (428 
nuc) following the standard sequencing procedures (Smip 
et al. 2013). We added to our dataset and aligned using Bio- 
Edit (HALL 1999) 251 cytochrome b sequences available in 
GenBank (Supplementary document I) to increase the tax- 
onomic and geographic scope. Our final phylogeny covered 
most of the geographic range of the Sand lizards except for 
the Caucasian L. a. ioriensis and the recently described 
L. a. mzymtensis. However, both forms are genetically in- 
distinguishable from L. a. boemica and L. a. exigua respec- 
tively (DORONINA 2021). To date the divergences within the 
L. agilis lineage we built a second dataset by assembling se- 
quences of mitochondrial 12s RNA and cytochrome b from 
Gene Bank (Supplementary document I) of lacertid lizards 
including all the representatives of the genus Lacerta and 
several outgroups. We included three combined sequences 
of the eastern (L. a. agilis) and western (L. a. exigua) line- 
ages as well as one L. a. garzoni. Then, we used some dated 
divergences between the outgroup species as points of cali- 
bration to date the divergences among these three L. agilis 
subspecies. The divergence between the Ereminae lacertids 
Mesalina and Acanthodactylus was placed at 16.0-13.1 Mya 
using a log-normal distribution (u = 2.7, © = 0.07, offset 
= o). Maximum divergence was established that matches 
the beginning of the aridification of the Saharo-Arabian 
region (FLOWER & KENNETT 1994), while the lower limit 
was established based on the estimated age of the first fos- 
sil record of Mesalina (RAGE 1976). Although this fossil was 
previously described as Eremias, this genus included at the 
time of the description Mesalina, geographically restrict- 
ed to Asia and temperate eastern Europe, making unlike- 
ly other assignation than this genus. The split between the 
two Iberian Peninsula Timon species was assumed to have 
occurred at 5.3-1.9 Mya based on their common divergence 
from the north African Timon species after the closure of 
the Strait of Gibraltar at the end of the Messinian salinity 
crisis (GAUTIER et al. 1994) and the oldest known fossil of 
Timon lepidus (Pleistocene, Acusti et al. 2010). This win- 
dow of divergence time was accounted for by a log-normal 
model setting u = 1.12, o = 0.24 with an offset = o. The last 
point of the calibration was the speciation event involving 
Lacerta trilineata and its sister lineage of green lizards that 
occurred 18.1-7.0 Mya, based on the first fossil record cer- 
tainly assigned to this species and this group (CERNANSKY 
& SYROMYATNIKOVA 2019). In this case, the log-model was 
fixed to u = 2.4, O = 0.23, offset=o. The 95% confidence in- 
terval of the three estimated divergences among sand lizard 
lineages was used to calibrate a phylogenetic tree of L. agilis 
based on cytochrome b sequences by setting log-normal 
distributions to account for the maximum and minimum 
values of the intervals. Bayesian phylogenetic analyses were 
carried out by Beast 2.6.2 (BOUCKAERT et al. 2019) under 
relaxed molecular clocks using the calibration schemes 
previously defined. We selected the GTR and HKY models 
with gamma for the analysis of cytochrome b plus 12sRNA 
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of Lacerta and outgroups after finding the best model for 
each gene and dataset utilizing jModeltest 2.1.10. (DAR- 
RIBA et al. 2012), and GTR for the cytochrome b of Sand 
lizards datasets. We ran three chains of 10° iterations tak- 
ing samples every 10‘ iterations for both analyses, check- 
ing for convergence and effective sample sizes using Tracer 
1.7.1. (RAMBAUT et al. 2018) and compiled the 100 last trees 
of each run to compute the maximum credibility tree in 
each analysis. We also analysed the genetic variability of 
each mitochondrial lineage by computing the haplotype 
and nucleotide diversities using DnaSP 6.12 (Rozas et al. 
2017). Biogeographic inferences were made by reconstruct- 
ing the evolution of ancestral ranges of the mitochondrial 
lineages of Sand lizards on the maximum credibility tree 
through BioGeoBears implemented in Rasp 4.2 (MATSKE 
2014). Thus, we compared the dispersal-extinction clado- 
genesis (DEC), Bayesian analysis of biogeography (BAYA- 
REALIKE) and the dispersal-vicariance (DIVA) models 
using Akaike Information Criterion corrected for small 
sample sizes (AICc, CAVANAUGH et al. 2019). The analy- 
sis was constrained to reconstruct ancestral areas formed 
by a maximum of two biogeographic units and those that 
are geographically adjacent, assigning each mitochondrial 
lineage these geographic units: Western Europe, Pyrenees, 
Carpathians, Balkans, South Crimea, Eastern Europe, and 
Caucasus. 


Species niche models 


We gathered data on the distribution of Lacerta agilis across 
its range from open databases and scientific papers (Sup- 
plementary document II). We also used our data from the 
Pyrenean populations based on fieldwork done by the au- 
thors (Amar et al. 2003). We evaluated the niche divergence 
in L. agilis of the main mitochondrial lineages found by the 
phylogenetic analysis. We calibrated statistical models using 
current relationships between lineage distribution and cli- 
mate conditions. These models were then used to infer niche 
suitability under fluctuating climatic scenarios, assuming 
that species responses are stable over time (MARTINEZ-MEy- 
ER et al. 2004, HOSNER et al. 2014, AMAT & ESCORIZA 2022). 

To build the best ecological niche models, we first gener- 
ated a smaller subset of the locations available in the data- 
bases, after removing those separated less than 10 km. This 
allows us to reduce spatial data autocorrelation and poten- 
tial model overfitting (YACKULIC et al. 2013). This process 
has been carried out using spThin routines (AIELLO-LAM- 
MENS et al. 2015) in R (R Core Development Team 2023). 
Niche models were built using Maxent 3.4.4 (PHILIPs et al. 
2006), a method suitable for generating projections based 
solely on presence data (WEsT et al. 2016). Several mod- 
els were built iteratively, testing successive combinations of 
feature classes (L, linear; Q, quadratic; H, hinge; P, prod- 
uct; T, threshold) and regularization multipliers and select- 
ing the optimal candidate using the modified Akaike cri- 
terion for finite sample sizes. The predictive capability of 
the models was evaluated using the area under the receiver 
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operating characteristic curve (AUC) and the continuous 
Boyce Index (CBI) (Dı Cota et al. 2017). CBI provides a 
better estimate of the predictive capacity of models based 
only on presence data; this index ranges between -1 and 1, 
with negative values indicating a wrong model and a val- 
ue close to 1 indicating almost perfect predictive accuracy 
(HIRZEL et al. 2006). These analyses were performed using 
ENMval functions (MUSCARELLA et al. 2014) in the R pro- 
gramming environment. 


Climate data 


We used 19 bioclimatic variables and a digital terrain mod- 
el available in WorldClim 2 (Fick & HIjMANS 2017) for 
niche modelling. Variable redundancy was estimated with 
the variable variance inflation factor (VIF) using a mod- 
el of increasing complexity (CRANEY & SURLES 2002). To 
calculate the VIF, we defined buffer polygons of 200 km 
around the species’ occurrences and generated 1000 ran- 
dom points within these polygons. The climate data ob- 
tained for the random points and the species’ occurrences 
were included in the logistic regression. We started with a 
minimal logistic model that included only BIO10 (mean 
temperature of the warmest quarter), BIOu1 (mean tem- 
perature of the coldest quarter), BIO12 (annual precipita- 
tion), and elevation. These variables were initially chosen 
due to their putative relevance in the life cycles of meso- 
philic temperate lizards, i.e. including summer and win- 
ter temperatures (egg embryonic development and hiber- 
nation) and a proxy for the amount of environmental hu- 
midity (VAN NULAND & STRIJBOSCH 1981, Li et al. 2013, Es- 
CORIZA & AMAT 2021). Elevation has also been included 
because the species’ southern populations are entirely con- 
fined to mountain habitats (NETTMANN et al. 1992, AMAT 
et al. 2003). Variables that showed a VIF > 10 have been ex- 
cluded from subsequent analyses (SCHROEDER et al. 1990). 

To assess niche suitability in the past, we selected the fol- 
lowing simulations: Pliocene, Piacenzian stage, (i) 3.3 Mya 
(cold) and (ii) 3.2 Mya (warm); Pleistocene, Calabrian stage 
(787 kya); Pleistocene, Chibanian stage (last interglacial, 
130 kya); Pleistocene, Last Glacial Maximum (21 kya); Pleis- 
tocene, Bolling—Allerod interstadial (14.7-12.9 kya); Pleis- 
tocene, Younger Dryas Stadial (12.9-11.7 kya); Holocene, 
Northgrippian (8.3-4.2 kya) (OTTO-BLIESNER et al. 2006, 
Dotan et al. 2015, HILL 2015, FORDHAM et al. 2017, BROWN 
et al. 2018, KARGER et al. 2021). These stages encompassed 
successive warm-cold oscillations, with peak-warm peri- 
ods (global temperatures 3°C higher than the current ones, 
DE LA VEGA et al. 2020) and cold-peak periods (global tem- 
peratures 5°C lower than current ones, KARGER et al. 2021). 
These projections were conducted using the whole set of the 
L. agilis species occurrences because realized niches could 
be limited by the presence of parapatric taxa and not only 
for their physiological responses (PANZACHI et al. 2015). 
Only in the case of the pair garzoni - agilis, we assessed 
their dynamics of niche evolution separately, because these 
contiguous taxa are completely allopatric. 
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Niche divergence and phylogenetic signal 


We also assessed the amount of niche divergence between 
pairs of L. agilis mitochondrial lineages. The relative niche 
position of each lineage was first visualized through prin- 
cipal component analysis (PCA) of the normalized envi- 
ronmental variables. This analysis was carried out in R (R 
Core Team 2023). 

Niche comparisons were conducted for all parapatric 
or contiguous lineages and allopatric and non-contigu- 
ous lineages only if they showed broad geographical dis- 
tributions. Niche divergence was assessed using several 
tests comparing the D index, which ranges between o (no 
overlap) to 1 (full overlap) (SCHOENER 1970). Range-break- 
ing tests determine sharp changes in environmental con- 
ditions at mitochondrial lineages’ boundaries. The linear 
test assumes that this limit between parapatric or contigu- 
ous lineages is linear, whereas the blob test assumes that it 
could be geometrically irregular, being more suitable when 
the lineage occupies geographical areas differing in extent 
(GLOR & WARREN 2011). Background tests assess whether 
the niches of the two species were more similar or dissimi- 
lar to each other than to the available conditions (WARREN 
et al. 2008). The statistical significance of these has been 


0.972 
0.634 j 
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obtained by comparing the observed values of taxa overlap 
with those simulated with 500 replicates, using the package 
ENMTools (WARREN et al. 2021). 

The test of phylogenetic effects over niche occupan- 
cy was evaluated using BLOMBERG’S K (BLOMBERG et 
al. 2003). BLOMBERG’S K values range between o and o, 
where the values of K < 1 indicate a less phylogenetic signal 
than that expected under Brownian motion (BLOMBERG et 
al. 2003). BLOMBERG’S K provides acceptable estimations 
even with small phylogenies, although it is susceptible to 
false positives (MUNKEMULLER et al. 2012). The value of 
K was calculated after 10,000 resamplings of the climate 
matrix, since it has been recommended not to use aver- 
age niche values in the estimation of the phylogenetic sig- 
nal (HARMON & Losos 2005). These calculations were con- 
ducted using the phytools package (REVELL 2012) in R. 


Results 


Our phylogenetic analysis placed sand lizards as a well- 
supported sister clade of that formed by Lacerta media, 
Lacerta pamphylica, and Lacerta trilineata (Fig. 1), estimat- 
ing the divergence between the two clades at 8.7-5.0 Mya. 
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Lacerta agilis garzoni 
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Lacerta media wolterstorffi 
Lacerta pamphylica 
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Lacerta trilineata cariensis 
Lacerta trilineata dobrogica 
Lacerta trilineata diplochondrodes 
Lacerta trilineata galatiensis 
Lacerta bilineata bilineata 
Lacerta bilineata chlorosecunda 
Lacerta bilineata chloronota 
Lacerta viridis meridionalis 
Lacerta viridis viridis 

Lacerta schreiberi 

Lacerta strigata 

Timon kurdistanicus 

Timon princeps 

Timon lepidus 

Timon nevadensis 

Timon pater 

Timon tangitanus 


Figure 1. Bayesian analysis of phylogenetic relationships of green lizards (Lacerta) based on a concatenated dataset of mitochondrial 
12s RNA and cytochrome b. Values on the nodes are the posterior probabilities. Divergences between taxa estimated using a relaxed 


molecular clock are shown in bars (95% confidence intervals). 
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Our results revealed that the origin of Sand lizards can be 
placed during the late Miocene-Pliocene within a time 
window of 3.8-1.8 Mya. Using the divergence dates among 
the three lineages of Sand lizards in this phylogenetic tree 
and the full dataset of cytochrome b sequences, we iden- 
tified several mitochondrial lineages with high support 
(Fig. 2). The most basal subclade (south Caucasian L. a. 
boemica) split at 2.9-1.7 Mya from the other L. agilis mito- 
chondrial subclades. The other group also split shortly af- 
ter (2.9-1.6 Mya): an eastern group constituted by the sub- 
species L. a. bosnica, L. a. exigua and L. a. tauridica, and a 
western group formed by L. a. agilis, L. a. chersonensis, L. a. 
argus, L. a. garzoni and an unnamed subclade confined to 
the Carpathian mountains. Within the eastern group, the 
Balkanic lineage diverged early (2.6-1.3 Mya) and became 
separated into three well-supported subclades, in the Pin- 
dos mountains, Bosnia, and the rest of the Balkans, with 
a minimum interclade divergence estimated between 
930-680 kya years ago. The other two lineages, the south 
Crimean L. a. tauridica, which was formed by two well- 
supported clades diverging at 1.2-0.4 Mya, and exigua, be- 
came mitochondrially distinct at 1.8-0.8 Mya. The western 
group began to diversify at 1.9-0.9 Mya originating the 
chersonensis lineage, and a group of lineages ranged from 


the easternmost, located in the Carpathians, to the west- 
ernmost, which included garzoni, agilis, and argus subspe- 
cies diverging at 1.5-0.8 Mya. Following our results, the 
origin of Pyrenean endemic garzoni was older than previ- 
ously stated, being at 1.2-0.7 Mya in its divergence from the 
common ancestor of the agilis plus argus subspecies, which 
diverged most recently (1.0-0.4 Mya). 

Sand lizard lineages differed in the amount of genetic di- 
versity (Table 1). Therefore, despite sampling thirty-five in- 
dividuals of the Pyrenean L. a. garzoni, this subspecies was 
uniform for cytochrome b in contrast to southern lineages, 
L. a. boemica, L. a. bosnica, and L. a. tauridica that showed 
higher haplotype and nucleotide diversity. The other line- 
ages had low levels of mitochondrial diversity even in the 
case of L. a. argus, L. a. agilis, and L. a. exigua with very ex- 
tensive geographic ranges. 

The reconstruction of the past biogeography of Lacerta 
agilis was based on the DIVA model, which performs bet- 
ter than the others (AICc weight = 0.670). We found the 
area formed by the Caucasus plus the adjacent part of east- 
ern Europe as the most likely ancestral range of Sand liz- 
ards. Evolving from the common ancestor, one lineage, L. a. 
boemica, remained confined to the eastern Caucasus, and 
another become separated into two groups and experienced 


Lacerta agilis agilis 


Lacerta agilis argus 


Lacerta agilis garzoni 


Carpathian clade 


Lacerta agilis chersoniensis 
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Figure 2. Bayesian analysis of cytochrome b sequences of 286 sand lizards. For simplicity, samples were collapsed into the nine main 
identified clades, most of them assignable to the currently recognised subspecies. The values of the nodes are the posterior prob- 
abilities. The bars represent the 95% confidence intervals of the divergence time among clades based on a relaxed molecular clock. 
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Table 1. Genetic variation of the main lineages of sand lizards 
based on mitochondrial cytochrome b. 


Lineage n Numberof Haplotype Nucleotide 
haplotypes diversity diversity 
argus 28 2 0.138+0.084 0.0006 
agilis 13 6 0.718+0.128 0.0055 
boemica 7 7 1.000+0.076 0.0082 
bosnica 19 10 0.825+0.084 0.0019 
carpathians 9 8 0.972+0.064 0.0085 
chersonensis 35 0.659+0.057 0.0036 
exigua 94 10 0.356+0.062 0.0023 
garzoni 36 1 0.000+0.000 0.0000 
tauridica 45 23 0.945+0.017 0.0071 


two independent step by step expansions to the rest of Eu- 
rope (Fig. 3). The westernmost group colonized indepen- 
dently the Pyrenees and Carpathians, from plain-dwelling 
ancestors, that also derived into the L. a. agilis, L. a. argus, 
and L. a. chersonensis subspecies. Another group colonized 
the Balkans and southern Crimea from a plain-dweller an- 
cestor that also evolved into L. a. exigua. The latter subspe- 
cies spread through the steppes towards Central Asia. 

The PCA showed the position of each lineage locality 
depending on climate variations (Fig. 4). The first axis (ex- 
plained variance = 41.02%) separated the lineages’ sites de- 
pending on the temperatures of the driest (BIOg9, variable 
partial contribution [vpc] = 23.9%) and the colder quarter 
(BIO11, vpc = 22.7%) and the amount of rainfall (BIO12, 
vpc = 15.7%) and its seasonality (BIO15, vpc = 15.2%) and 
the second axis (explained variance = 24.59%) represent- 
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ed the altitudinal gradient (altitude, vpc = 36.7%) and the 
temperatures of the warmest period (BIO10, vpc = 25.7%) 
(Fig. 4). 

The models performance was good, according to the 
AUC and CBI, and the mapped predictions were reliable 
(Supplementary document III). Under current conditions, 
the model forecasts showed broad regions climatically suit- 
able distributed sequentially from west to east, extending 
across the Great European Plain to the West Siberian Plain 
(agilis-argus-chersonensis-exigua) with smaller suitable 
regions restricted to the Pyrenean, Carpathian, Dinaric 
Alps, Pindus, Tauric (Crimean), and Caucasus mountain 
systems, which also sequentially correspond to the endem- 
ic montane subspecies (Fig. 5). 

Paleoclimatic models showed important fluctuations 
between the gain-loss of suitable habitats in Europe dur- 
ing the last 3.3 Mya. In the Pliocene, the mapped projec- 
tions revealed a significant loss of suitability across the 
Great European Plain and the Mediterranean peninsulas, 
when at this moment the suitable regions were confined 
to the great mountainous systems and the north-eastern 
margins of the Great European Plain (Fig. 6). During cold- 
er interstadial periods, the Great European-West Siberian 
Plains become suitable for this species, but not most of the 
Mediterranean peninsulas. During extremely cold periods, 
such as the last glacial maximum (21 kya), there was a ma- 
jor loss of suitable habitats in the northern and eastern dis- 
tribution, including the West Siberian Plain and most of 
the Great European Plain, while the Atlantic margins of the 
Great European Plain and the north of the Mediterranean 
Peninsulas likely act as refugia (Fig. 6). Interestingly, habi- 
tats at the glacier margins could have been suitable for this 
species in western Europe. 


Figure 3. Reconstruction of the historical biogeography of sand lizards based on the dispersal-vicariance model. The image depicts 
the inferred process of dispersal to western Europe experienced by the lineage comprising the clades assignable to the garzoni, agilis, 
argos, chersonensis subspecies and the unnamed clade inhabiting the Carpathians (in blue) and the formed by bosnica, tauridica, and 
exigua (in red). The Caucasus and adjacent plains of eastern Europe probably were the ancestral range of Lacerta agilis. The boemica 
subespecies (in black) was the first clade to evolve independently to the others and lives in the Eastern Caucasus. 
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PC1 (41.02%) 


Figure 4. PCA scatter plot representing the relative position of each lineage of Lacerta agilis depending on climatic variation. The 
localities of each species are shown encompassed by a convex hull. ALT, terrain elevation. BIO, climatic variables. 


dz — ss: D 


Figure 5. Ecological niche models built with Maxent based on WorldClim 2 variables representing the niche of Lacerta agilis lineages, 
generated from Worldclim bioclimatic variables and terrain elevation. (A) Lacerta agilis agilis; (B) L. a. argus; (C) L. a. garzoni; (D) L. a. 
“Carpathians”; (E) L. a. bosniensis; (F) L. a. chersonensis; (G) L. a. exigua; (H) L. a. tauridica; (1) L. a. boemica. 
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In the case of L. a. garzoni, the Pyrenean had suitable 
habitats for this subspecies since 787 kya. Pleistocene mod- 
els indicated the persistence of suitable climatic conditions 
for L. a. garzoni on the southeastern edges of the Pyrenees, 
maintaining geographic isolation from the nominate sub- 
species (Fig. 7). 

The lack of statistical significance for the comparisons 
of the linear and blob range-breaking tests allows us to re- 
ject the hypothesis of sharp environmental variations at the 
subspecies range boundaries (Table 2). The background 
tests revealed higher similarity than expected by chance in 
L. a. agilis - L. a. argus and L. a. bosnica - L. a. ‘Carpathi- 
ans, and lower similarity than expected by chance in L. a. 
agilis - L. a. garzoni, L. a. agilis - L. a. exigua, L. a. argus - 
L. a. ‘carpathians; and L. a. argus and the two eastern plains 
subspecies (L. a. chersonensis and L. a. exigua) (Table 3). 

The phylogenetic signal (BLOMBERG’s K) varied wide- 
ly depending on the environmental variable evaluated, 
but tended to be highest (i.e. more conserved variables) 
in those describing extreme temperatures: warmest (mean 
K = 0.825) and coldest quarter (mean K = 0.799) and pre- 


-40 
12.9 kya 


-10 o 


Table 2. Ecological niche hypothesis tests: p-values of the range 
break tests linear and blob reject. Non-significant values rejected 
the hypothesis of sharp environmental variations at the species 
range boundaries. AL, allopatric and non-contiguous subspecies. 


Linear Blob 
L. a. agilis - L. a. argus 0.373 0.078 
L. a. agilis - L. a. garzoni 0.412 0.129 
L. a. agilis - L. a. chersonensis AL AL 
L. a. agilis - L. a. exigua AL AL 
L. a. argus - L. a. bosnica 0.392 0.231 
L. a. argus - L. a. ‘carpathians’ 0.019 0.255 
L. a. bosnica - L. a. ‘carpathians’ 0.333 0.178 
L. a. argus - L. a. chersonensis 0.137 0.337 
L. a. argus - L. a. L. a. exigua AL AL 
L. a. ‘carpathians’ - L. a. chersonensis 0.196 0.406 
L. a. exigua - L. a. chersonensis 0.159 0.149 
L. a. exigua - L. a. boemica 0.294 0.099 
L. a. exigua - L. a. tauridica 0.082 0.079 


` 7 
-10 0 3 50 6 


Figure 6. Paleo-climatic projections, based on models statistically calibrated with the current climatic conditions: Pliocene, Piacenzian 
stage, (i) 3.3 ma (cold) and (ii) 3.2 ma (warm); Pleistocene, Calabrian stage (787 kya); Pleistocene, Chibanian stage (last interglacial, 
130 kya); Pleistocene, Last Glacial Maximum (21 kya); Pleistocene, Bølling-Allerød interstadial (14.7-12.9 kya); Pleistocene, Younger 
Dryas Stadial (12.9-11.7 kya); Holocene, Northgrippian (8.3-4.2 kya). 
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Table 3. Ecological niche hypothesis tests: Schoener’s D and 
p-values of the background tests. Non-significant (n.s.) values 
rejected the hypothesis of greater similarity or divergence than 


would be expected by chance. 


Direction 
L. a. agilis - Empirical 0.225 Similarity 
L. a. argus Mean simulated 0.165 
p 0.020 
L. a. agilis - Empirical 0.015 Dissimilarity 
L. a. garzoni Mean simulated 0.248 
p 0.019 
L. a. agilis - Empirical 0.163 n.s 
L. a. chersonensis Mean simulated 0.185 
p 0.078 
L. a. agilis - Empirical 0.234 Dissimilarity 
L. a. exigua Mean simulated 0.254 
p 0.039 
L. a. argus - Empirical 0.208 n.s 
L. a. bosnica Mean simulated 0.201 
p 0.451 
L. a. argus - Empirical 0.363 Dissimilarity 
L. a. ‘carpathians Mean simulated 0.723 
p 0.019 
L. a. bosnica - Empirical 0.450 Similarity 
L. a. ‘carpathians’ Mean simulated 0.219 
p 0.019 
L. a. argus - Empirical 0.276 Dissimilarity 
L. a. chersonensis Mean simulated 0.629 
p 0.048 
L. a. argus - Empirical 0.183 Dissimilarity 
L. a. exigua Mean simulated 0.392 
p 0.019 
L. a. ‘carpathians - Empirical 0.083 n.s 
L. a. chersonensis Mean simulated 0.128 
p 0.137 
L. a. exigua - Empirical 0.204 n.s. 
L. a. chersonensis Mean simulated 0.224 
p 0.157 
L. a. exigua - Empirical 0.209 Dissimilarity 
L. a. boemica Mean simulated 0.279 
p 0.019 
L. a. exigua - Empirical 0.067 Dissimilarity 
L. a. tauridica Mean simulated 0.177 
p 0.039 


cipitation seasonality (mean K = 0.792) 


. However, we 


found a lower phylogenetic signal (K < 0.707) in the oth- 
er variables of the environmental niche (terrain elevation, 
annual mean precipitations and mean temperatures of the 
driest and wettest periods of the year) (Table 4). 


Discussion 
Our study infers a realistic biogeographic scenario to ex- 
plain how Sand lizards achieved their current huge geo- 


graphic range, emphasizing the benefits of combining phy- 
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Table 4. Estimation of the phylogenetic signal (BLoMBERG’s K) 
in the diversification of the ecological niche in Sand lizards. The 
mean value and the standard error after 10,000 resamplings of 
the possible values of the niche for each subspecies are shown. 


Code Name BLOMBERG'’s K 
Elevation 0.707+0.002 
BIO8 Temperature wettest quarter 0.671+0.001 
BIO9 Temperature driest quarter 0.706+0.001 
BIO10 ‘Temperature warmest quarter 0.825+0.002 
BIO11 ‘Temperature coldest quarter 0.799+0.001 
BIO12 Annual precipitation 0.685+0.001 
BIO15 Precipitation seasonality 0.792+0.002 


logeographic and paleoclimatic methods (e.g. WIELSTRA et 
al. 2013). 

In general, as expected, climate models indicated that 
the geographical range of the whole species has not been 
stable in the last 3 Mya. The models also suggested that 
extensive secondary contacts between the plain-dwelling 
subspecies have only occurred during the interglacial-in- 
terstadial periods. The absence of abrupt discontinuities in 
the limits of the subspecies and niche similarity between 
several pairs of geographically contiguous subspecies sug- 
gests allopatric differentiation processes, compatible with 
the use of separated refugia during prolonged climatically 
hostile periods. 

Our results revealed that this species has spread through 
much of the cool-temperate regions of Eurasia favoured by 
a certain plasticity of the climatic niche, particularly in the 
amount and seasonality of precipitations. This plasticity 
has allowed these subspecies to colonize habitats exposed 
to very seasonal hydrological regimes, such as the cold 
steppes of Central Asia, or not, such as the oceanic climates 
of the European Atlantic plain (ANANJEVA et al. 2006). On 
the contrary, the tolerance to temperatures (particularly 
to the maximum ones) is more conserved, which would 
have limited the penetration of this species in regions with 
a Mediterranean climate, despite the numerous potential 
ways of dispersal available during their long evolutionary 
history. These findings mirror those of previous studies in 
the Mediterranean region, including several species and 
genera, suggesting that maximum temperatures are an im- 
portant factor that limits the ecological diversification in 
lacertid lizards (Escoriza et al. 2021). 

The similarity of the niches when comparing closely re- 
lated subspecies and the absence of sharp environmental 
discontinuities, such as L. a. agilis and L. a. argus and L. a. 
chersonensis and L. a. exigua, suggests that range bounda- 
ries are likely to be constrained by competitive interactions 
between conspecific populations. The results showed that 
niche divergence is high between some subclades with deep 
genetic divergence, i.e when L. a. boemica is compared to 
the surrounding plain-dweller subspecies L. a. exigua. In 
this case, niche divergence may favour their separation 
even if range-breaking tests did not support abrupt changes 
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in their range boundaries. However, in this case, subspecies 
replacement can occur at a lower spatial resolution, associ- 
ated with gradual, fine-scale changes in the landscape (e.g. 
riparian forest to steppe), not captured in our study. 

Based on our results, it is likely that the lineage of 
L. agilis arose in the region comprising the Caucasus and 
the adjacent plains between the Pleistocene—Pliocene tran- 


L. agilis garzoni 
21 kya 


sition. Therefore, the origin of the species may be more 
recent than previously stated (Pliocene, YABLOKOV et al. 
1980, KALYABINA et al. 2001, JOGER et al. 2006) and is not 
clear whether it evolved from a plain-dweller or mountain 
ancestor. In this sense, the basal position of the boemica 
subspecies respect to the other lineages could equally indi- 
cate a Caucasian origin of sand lizards or an early coloniza- 


TE 


L. agilis agilis 
" 21 kya 


Figure 7. Detailed paleo-climatic projections of the Pyrenees, based on models statistically calibrated with the current climatic condi- 
tions specifically for L. a. agilis and L. a. garzoni: Pleistocene, Last Glacial Maximum (21 kya); Pleistocene, Bølling-Allerød interstadial 
(14.7-12.9 kya); Pleistocene, Younger Dryas Stadial (12.9-11.7 kya); Holocene, Northgrippian (8.3-4.2 kya). 
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tion of this mountain range. The split between boemica and 
the common ancestor of the western and eastern lineages 
is old enough to make this subspecies the candidate to be 
recognised as a full species (JOGER et al. 2006). However, 
the taxonomic revision of Sand lizards will be greatly ben- 
efited from an analysis of nuclear DNA markers to quan- 
tify the amount of introgression between boemica and the 
other lineages. 

Paleoclimatic projections combined with molecular di- 
vergence times rejected previous hypotheses that proposed 
a Holocenic origin of the western lineage (KALYABINA et 
al. 2001). In this sense, paleoclimatic projections indicat- 
ed extensive habitat suitability across the plains of western 
Europe at a minimum of 787 kya. This period encompass- 
es the lower bounds of estimated divergence times of the 
agilis, garzoni, and argus subspecies. Our results rejected a 
recent colonization of the Pyrenees (PALACIOS & CASTRO- 
VIEJO 1975, CARRETERO et al. 2002), reinforcing the sub- 
specific identity of the populations that live in the massif. 

The mapped projection showed larger habitat suitabil- 
ity for the garzoni subspecies, than the actual distribution 
of the species. This incongruence could be caused by sev- 
eral reasons: the use of a very small glacial refugium, which 
would not have allowed the re-expansion along the rap- 
idly emerging favourable habitats, the presence of impass- 
able barriers, or competition with other congeneric lacer- 
tid species (HELTAI et al. 2015, RYAN & GUNDERSON 2021). 
Paleoclimatic projections indicated that during the Last 
Glacial Maximum, suitable conditions for this subspecies 
were restricted to a very small region on the southeast- 
ern edge of the Pyrenees. A putative glacial refuge could 
be located in the valley of the Tec River, thus facilitating 
the connection of the two main groups of Pyrenean pop- 
ulations (ARRIBAS 1999, AMAT & ROIG 2003, GENIEZ & 
CHEYLAN 2012). This valley is the only one with an east- 
east orientation inhabited by the species in the Pyrenees, 
and this orientation (toward the sea) may favour milder 
thermal conditions. The existence of only one glacial refu- 
gium for Pyrenean populations is congruent with the ab- 
sence of mitochondrial DNA variation in our sampling, 
similar to other microendemic lizards (Iberolacerta aure- 
lioi and I. aranica, CARRANZA et al. 2004). Other Pyrene- 
an species, more widely distributed, however, have a more 
complex geographic structure (Zootoca vivipara, MILA et 
al. 2013 and Iberolacerta bonnali, CARRANZA et al. 2004, 
FERCHAUD et al. 2015), suggesting multiple climate refu- 
gia through the Pyrenees. The Last Glacial Maximum pro- 
jection also predicted large habitat suitability for the agilis 
subspecies along the periphery of the main chain that over- 
lapped with the estimated for garzoni. However, the habi- 
tat suitability of the intermediate region between the Pyr- 
enees and the Central Massif was low, suggesting a main- 
tained isolation between the populations of both mountain 
ranges. Therefore, the absence of this species in the central 
and western Pyrenees could be explained by a slow process 
of expansion from this refugium (or others located on the 
southeastern edges of the Pyrenees). Although competi- 
tion with Lacerta bilineata could also be arguable, Sand liz- 
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ards roughly coexist with other green lizards (Lacerta viri- 
dis, HELTAI et al. 2015). 

Our mapped projections indicated that the great plains 
of Europe become corridors for an early dispersal, likely 
from the plains close to the Black Sea to colonize indepen- 
dently Balkan mountain ranges, than was previously hy- 
pothesised, but also rejecting the role of the Balkans as the 
centre of dispersal for the eastern lineage (KALYABINA et 
al. 2001). The paleoclimatic projections indicated that the 
expansion towards Central Asian steppes by L. a. exigua 
occurred recently, around 8,300 years (also suggested by 
KALYABINA et al. 2001). 

Noticeably, we found higher mitochondrial diversity in 
all southern montane clades with the sole exception of the 
Pyrenean garzoni due to the prolonged stable habitat suita- 
bility in the great mountain ranges of the Balkans and east- 
ern Europe. In contrast, paleoclimatic projections and low- 
er mitochondrial diversity suggested a smaller extent of the 
refugia of exigua, chersonensis, argus, and agilis subspecies 
with few populations. This scenario was reversed, after the 
glacier retreat allowed the expansion to the plain-dweller 
lineages across the Eurasian plains. 

Our study evidenced that Pleistocene glacial cycles have 
deeply shaped the evolutionary history of L. agilis, leading 
to the isolation of mitochondrial lineages in the southern 
Europe refugia, in plains, and separated mountain rang- 
es. Subsequent range expansions could be more rapid in 
the plain-dweller subspecies (agilis, argus, chersonensis, 
and exigua), enabling them to expand over vast geographic 
areas, while mountain lineages remained more confined. 
Interestingly, except for the Pyrenean subspecies, south- 
ern montane lineages have a higher mitochondrial diver- 
sity more broadly distributed than plain-dwelling lineages, 
suggesting larger glacial refugia in the former or, at least, 
that the properties of the habitats allowed larger popula- 
tion sizes. 
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